In recent years, the research of fractional calculus has exercised strong attractions.^[1–10] In the field of engineering, a lot of problems can be better described by fractional-order equations, for instance, the description of the memory and hereditary properties in various materials and processes,^[9,10] the model of rheological medium and viscoelastic fluid.^[11–13] Many works had been fulfilled on fractional-order dynamic system, where the dynamic properties and the effects of the fractional-order terms are concerned. Those classical quantitative analysis methods are extended successfully to analyze the dynamical properties of fractional-order systems. Wahi and Chatterjee^[14] demonstrated the method of averaging for conservative oscillators which may be strongly nonlinear, under small perturbations including delayed and/or fractional derivative terms. Chen et al.^[15] derived the averaged Itô equations of amplitude modulation and phase difference via stochastic averaging method. Shen et al.^[16] investigated the primary resonance of dry-friction oscillator with fractional-order PID controller of velocity feedback by Krylov–Bogolubov–Mitropolsky (KBM) asymptotic method. Xu et al.^[17] investigated stochastic Duffing oscillator with fractional-order damping term by combining with Lindstedt–Poincaré (LP) method and the multiple-scale approach. Shen et al.^[18] extended the incremental harmonic balance (IHB) method to analyze the dynamical properties of fractional-order nonlinear oscillator, which is a very effective method to deal with both strongly and weakly nonlinear systems. Xie and Lin^[19] investigated the free vibration of van der Pol oscillator with small fractional damping by using the method of two-scale expansion. Rand et al.^[20] obtained the approximation of fractional Mathieu equation by the harmonic balance method. Yang and Zhu^[21] obtained the approximate solutions of the fractional-order linear system with harmonic excitation by the harmonic balance method. Guo et al.^[22] proposed an analytical technique based on the method of harmonic balance for predicting and generating the steady-state solution of the fractional differential system.

The Duffing–van der Pol equation is essentially a combination of Duffing and van der Pol equations, and it can be used as a mathematical model in many fields,^[23] such as aircraft wings at high angles of attack, thin panels in supersonic flows, and the structure vibration caused by the flow, etc. Lots of research work has been carried out on the Duffing–van der Pol system. Awrejcewicz and Mrozowski^[24] studied chaotic dynamics of a particular non-linear oscillator having Duffing type stiffness, van der Pol damping and dry friction. Kimiaeifar et al.^[25] solved the Duffing–van der Pol equation analytically using homotopy analysis method. Kyziołand Okniński^[26] obtained the periodic steady-state solutions of the Duffing–van der Pol equation by the KBM approach. Recently, Duffing–van der Pol oscillator with fractional-order derivative had attracted more and more attentions. For example, Leung et al.^[27] investigated a Duffing–van der Pol oscillator having fractional derivatives and time delays by the residue harmonic method. Matouk^[28] studied the stability analysis of the fractional-order modified autonomous van der Pol–Duffing circuit by using the fractional Routh–Hurwitz criteria.

The harmonic balance (HB) method is a widely used method to deal with both integer order systems and fractional-order systems with strong and weak nonlinearity.^{[20,21,29–31]} Accordingly, the Duffing–van der Pol oscillator with fractional-order damping term is studied by HB method in this paper. The paper is organized as follows. Section 2 presents the amplitude–frequency response equation of primary resonance by using HB method. In Section 3, the stability condition for steady-state solution is presented. In Section 4, the numerical solutions of system are obtained by the power series method. In Section 5, the bifurcations of primary resonance are investigated. The main conclusions are drawn at last.

The fractional-order Duffing–van der Pol oscillator considered is as follows: 1 where m is the mass, c₁ and c₂ the fractional-order damping coefficients, k is the linear stiffness coefficient, α₁ is the nonlinear stiffness coefficient, F₁ is the amplitude of the external harmonic excitation, ω is the angular frequency of the external harmonic excitation. And p is the order of fractional-order derivative, which is restricted as 0 ≤ p ≤ 1. The Caputo’s definition of fractional-order derivative is adopted here. The fractional-order term in Eq. (1) is defined as 2 where Γ (z) is Gamma function satisfying Γ (z + 1) = zΓ (z).

To study the approximate periodic solutions by HB method, and equation (1) can be transformed into 3 In Eq. (3) x″ denotes the second derivative with respect to t, μ₁ = c₁/m, μ₂ = c₂/m, , α = α₁/m, F = F₁/m.

The first-order approximate periodic solution contains only one harmonic term, which could be supposed as 4 where A is the amplitude and θ is the phase to be determined.

By introducing φ = ωt + θ, one could obtain a transformation of Eq. (3) 5 where denotes the second derivative with respect to φ.

Substituting Eq. (4) into the left-hand side of Eq. (5), for the integer-order parts one could obtain 6 For the fractional-order part in Eq. (5), if p = 1 one could have 7

When 0 ≤ p < 1, the fractional-order term in Eq. (5) can be taken as a periodic function with the period tends to infinity. Expanding into Fourier series and taking N-th harmonic terms, one could have 8a 8b 8c 8d

In order to calculate the Fourier coefficients, two important results about the definite integrals in Ref. [18] are introduced as 9a 9b

Firstly, E_1r is calculated. Substituting s = φ − τ, ds = −dτ into Eq. (8c), one could obtain 10a 10b 10c

Using integration by parts and Eq. (8c), one could obtain 11a 11b 11c

Substituting Eq. (11) into Eq. (10), one could have 12

Similarly, one could obtain 13a 13b

Furthermore, one could obtain 14

From Eq. (7) and Eq. (14), it can be concluded that the fractional-order part in Eq. (5) can be expressed as the unified form by Eq. (14) for the fractional order 0 ≤ p ≤ 1.

Substituting Eq. (6) and Eq. (14) into Eq. (5) and remaining first-order harmonic terms, one could have 15

According to harmonic balance principle, one could obtain 16a 16b

Eliminating the parameter θ and substituting the original parameters of system, the amplitude–frequency response equation can be obtained as 17 where And the phase-frequency response equation can be expressed as 18

From Eq. (17), the peak value of the primary resonance can be obtained as 19 when K = mω².

In Eq. (17) and Eq. (18), the parameter K denotes the equivalent stiffness coefficient of the system. The parameter K(p) is defined as the equivalent stiffness coefficient, and C(p) is defined as the equivalent damping coefficient, which are both induced by the fractional-order term. Letting , and ω_r can be defined as the equivalent natural frequency of system. When p = 1, one could get K(p) = 0 and C(p) = (c₂A²/4) − c₁, it is obviously consistent with the integer–order case. It can be found that the viscoelastic damping has equivalent stiffness and damping effect at the same time.

The stability of the steady-state solution is studied in this section. The steady-state solution can be expressed as 20 And x_p satisfies Eq. (16). The small perturbations near the periodic motion x_p can be assumed as ΔA and Δθ. Let and . Substituting them into Eq. (16) and ignoring higher-order terms, it yields 21a 21b Substituting the original parameters of system, one could have 22a 22b Eliminating from Eq. (21) based on Eq. (16), and the characteristic determinant can be obtained as 23 where A₁ = −C(p)ω, , ,

Expanding the determinant, the characteristic equation can be obtained as 24

From Eq. (23), the stability condition for the steady-state solution can be expressed as 25 It can be concluded that the solution in the middle is unstable, when the amplitude–frequency response curve has multiple solutions.

In order to verify the precision of the approximate analytical solution, the numerical results are also provided. The relationship about the explicit numerical approximation of the power series method in Ref. [2] is 26 where t_n = nh is the sample points, h is the sample step, and is the fractional binomial coefficient with the iterative relationship as 27 According to Eq. (26) and Eq. (27), the numerical scheme for Eq. (1) can be expressed as 28a 28b 28c The sample step is set as h = 0.004, and the total computation time is 200 s. After omitting the temporary response in frontal 120 s, the peak value of the posterior response is taken as the steady-state amplitude by the numerical results.

In order to investigate the dynamic response characteristics of the system, the demonstration system parameters are defined as system parameters are defined as m = 10, k = 40, c₁ = 0.2, c₂ = 0.5, α₁ = 0.5, p = 0.5, F₁ = 2. Based on Eq. (17), one could analytically obtain the amplitude–frequency response curve, which is shown in Fig. 1 by the solid points. According to Eq. (28), one could obtain the amplitude–frequency response by numerical integration, which are also shown in Fig. 1 and denoted by the circles. From the observation of Fig. 1, it could be found that the approximate analytical solution agrees well with the numerical results.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Comparison of amplitude–frequency response.

In this section, the bifurcations of primary resonance for the given parameters and ω = 2.1 are presented. When the parameters of the system are changed respectively, the behaviors of resonant amplitudes can be obtained by analytical solution based on Eq. (17), and the results are shown in Figs. 2–6 by the solid points. From these figures, it can be seen that there exist multi-solutions and jump phenomenon all of them, and the intermediate solution is an unstable solution. For the parameters F₁, α₁, p, and c₂, the resonant amplitude maybe has only single solution by choosing appropriate parameter values.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Resonant amplitude with F1 for m = 10, k = 40, c1 = 0.2, c2 = 0.5, α1 = 0.5, p = 0.5, ω = 2.1.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Resonant amplitude with α1 for m = 10, k = 40, c1 = 0.2, c2 = 0.5, F1 = 2, p = 0.5, ω = 2.1.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Resonant amplitude with p for m = 10, k = 40, c1 = 0.2, c2 = 0.5, α1 = 0.5, F1 = 2, ω = 2.1.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Resonant amplitude with c1 for m = 10, k = 40, c2 = 0.5, α1 = 0.5, F1 = 2, p = 0.5, ω = 2.1.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Resonant amplitude with c2 for m = 10, k = 40, c1 = 0.2, α1 = 0.5, F1 = 2, p = 0.5, α = 2.1.

From Eq. (19), it can be concluded that the peak value of the primary resonance is related with the amplitude of the external harmonic excitation F₁, the equivalent natural frequency ω_r, the order of fractional-order derivative p. The peak value of the primary resonance trends to increasing along with the F₁ increased and decreased. When the order of fractional-order derivative is set as p = 1, equation (1) becomes a traditional integer–order system. In order to make the comparison more obvious, the amplitude of the external harmonic excitation is set as F₁ = 10. The amplitude–frequency response curves for fractional-order and integer–order systems can be obtained according to Eq. (17), which are shown in Fig. 7. It can be seen that, the peak resonant amplitude of the integer–order system is smaller than that of the fractional-order system.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Amplitude–frequency response curves for m = 10, k = 40, c1 = 0.2, c2 = 0.5, α1 = 0.5, F1 = 10.

The primary resonance of fractional-order of Duffing–van der Pol oscillator is studied in this paper. The approximate analytical solution is obtained by HB method, and the stability condition for steady-state solution is presented based on Lyapunov theory. By comparing with the numerical solutions, the accuracy of the analytical results is demonstrated. Based on amplitude–frequency response equation, the bifurcations of primary resonance are analyzed, and the effect on amplitude–frequency characteristic for the order of fractional-order term is pointed out. It may be helpful to analyze the forced vibration of other similar fractional-order oscillators.